第十六讲 R-双向方差分析1

在“R与生物统计专题”中，我们会从介绍R的基本知识展开到生物统计原理及其在R中的实现。以从浅入深，层层递进的形式在投必得医学公众号更新。

前面第十四和十五讲，我们介绍了单因素方差分析及其检验假设的验证（第十四讲 R-单因素方差分析1、第十五讲 R-单因素方差分析2）。

当我们针对一个因素进行分组时，单因素方差分析可以比较多组间有无差异。但是，当我们存在随机区组设计实验（又称配伍组设计）的时候，单因素方差分析便无法满足要求。于是需要用到今天介绍的：双向方差分析。

这将分两小讲进行解说，第十六讲主要解说双向方差分析的概念和前期数据展现及假设条件验证；第十七讲将介绍具体针对平衡设计的实验和不平衡设计的实验的方差分析的R实现。

随机区组设计

首先，我们先来熟悉一个概念，随机区组设计（randomized block design），又称配伍组设计，它是配对设计的扩展。

具体做法是：先按照实验结果的非处理因素（如性别、体重、年龄等）将实验对象配对成区组（block），再分别将各区组的实验对象随机分配到处理组或者对照组（引自《医学统计学》第四版，人民卫生出版社）。

比如，将小鼠按体重匹配的原则，随机分成5个区组，每个区组20只小鼠。再将20只小鼠按2:2随机分配到处理组和对照组。在这个过程中，我们涉及了两个因素，一个是按体重分配到区组，一个是处理与否。

1. 双向方差分析概念

双向方差分析（two-way analysis of variance test，ANOVA）用于同时评估两个分组变量（A和B）对响应变量的影响。

分组变量

分组变量也称为因素（factor），比如区组。

一个因素的不同类别（组）称为水平（level），比如处理组和不处理两个水平。

水平数可以在各个因素之间变化。因素的水平组合称为单元（cell），即每个区组中的每个处理／对照组。

• 当单元内的样本大小相等时，我们具有所谓的平衡设计。

  在这种情况下，可以应用标准的双向方差分析。

  比如，5个区组中每个区组内的处理组小鼠和对照组中小鼠都是10只。

• 当自变量每个级别内的样本量不相同时（不平衡设计的情况），应采用不同的方差分析方法。

  比如，区组间处理组和对照组小鼠数量不一样多。

2. 双向方差分析检验假设

1. 因素A的均值没有差异

2. 因素B的均值没有差异

3. 因素A和B之间没有相互作用

情况1和2的备择假设是：均值不相等。

情况3的备择假设是：A和B之间存在相互作用。

3. 双向方差分析的假设

与所有方差分析检验一样，双向方差分析假设每个单元内的观测值呈正态分布并且满足方差齐性。这一部分请参看往期推送请参看第十五讲 R-单因素方差分析2。

4. 双向方差分析：平衡设计

平衡设计即独立分组水平内的样本数量相等的情况。比如，5个区组中每个区组内的处理组小鼠和对照组中小鼠都是10只。

4.1将数据导入R

在这里，我们将使用名为ToothGrowth的内置R数据集。它包含一项评估维生素C对豚鼠牙齿生长的影响的研究数据。实验在60只豚鼠上进行，其中每只豚鼠通过两种递送方法（橙汁，OJ，或抗坏血酸，VC）分别接受三种剂量水平的维生素C量（0.5、1和2 mg /天， VC）。实验者测量了牙齿生长的长度，数据示例如下所示。

# 导入R内自带的ToothGrowth数据集
library(datasets)
data(ToothGrowth)

# 将数据存储在变量my_data中
my_data <- ToothGrowth

4.2 检查你的数据

# 显示前六行
head(my_data)

输出结果

len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

# 显示数据的结构
str(my_data)

输出结果

'data.frame':   60 obs. of  3 variables:
$ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
$ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
$ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

从上面的输出中，R将“剂量（dose）”视为数值型变量。我们将其转换为因素变量（即分组变量），如下所示。

# 将因素转换为水平变量
# 将剂量0.5，1，2分别转换为 "D0.5", "D1", "D2"
my_data$dose <- factor(my_data$dose,
levels = c(0.5, 1, 2),
labels = c("D0.5", "D1", "D2"))
head(my_data)

输出结果

len supp dose
1  4.2   VC D0.5
2 11.5   VC D0.5
3  7.3   VC D0.5
4  5.8   VC D0.5
5  6.4   VC D0.5
6 10.0   VC D0.5

研究问题：牙齿的长度是否取决于供应方式和剂量。

• 生成频率表：

table(my_data$supp, my_data$dose)

输出结果

D0.5 D1 D2
OJ   10 10 10
VC   10 10 10

这是一个2X3设计单元，其中因素是供应方式supp和剂量dose，每个单元中有10个研究对象。这是一个平衡的设计。

4.3可视化您的数据

箱形图和折线图可用于可视化组间差异：

• 箱形图用于按两个因素的水平对数据进行分组。

• 双向交互作用图，它绘制因素的双向组合的响应的平均值（或其他内容），从而说明可能的交互作用。

• 我们将使用ggpubr R包进行基于ggplot2的简单数据可视化。

• ggpubr的安装请参看：第五讲 R-数据描述性统计分析作图

4.3.1 箱形图

library("ggpubr")
ggboxplot(my_data, x = "dose", y = "len", color = "supp",
palette = c("red", "blue"))

4.3.2 双向交互作用图

library("ggpubr")
ggline(my_data, x = "dose", y = "len", color = "supp",
add = c("mean_se", "dotplot"), #添加均数和标准差
palette = c("red", "blue"))

4.3.3 其他绘图方式：

# 箱形图
boxplot(len ~ supp * dose, data=my_data, frame = FALSE,
col = c("red", "blue"), ylab="Tooth Length")

# 双向交互作用图
interaction.plot(x.factor = my_data$dose, trace.factor = my_data$supp,
response = my_data$len, fun = mean,
type = "b", legend = TRUE,
xlab = "Dose", ylab="Tooth Length",
pch=c(1,19), col = c("red", "blue"))

用于函数interact.plot()的参数：

• x.factor：

  要在x轴上绘制的因素：

  dose。

• trace.factor：

  要绘制为线的因素：

  sup

• response：

  给出响应的数值变量：

  len

• type：

  图的类型。

  允许的值包括p（仅用于点），l（仅用于线）和b（用于点和线）。

好了，本期讲解就先到这里。小伙伴们赶紧试起来吧。

在之后的更新中，我们会进一步为您介绍R的入门，以及常用生物统计方法和R实现。欢迎关注，投必得医学手把手带您走入R和生物统计的世界。

提前打个预告，在下一讲中，我们将为您详细介绍双向方差分析在平衡设计和不平衡设计的配对区组实验中的具体R语言实现。

当然啦，R语言的掌握是在长期训练中慢慢积累的。一个人学习太累，不妨加入“R与统计交流群”，和数百位硕博一起学习。

快扫二维码撩客服，

带你进入投必得医学交流群，

让我们共同进步！

↓↓

- END -

长按二维码关注「投必得医学」，更多科研干货在等你！